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A thorough study of shear stress within the lattice Boltzmann method is provided. Via standard 
multiscale Chapman-Enskog expansion we investigate the dependence of the error in shear stress on 
grid resolution showing that the shear stress obtained by the lattice Boltzmann method is second 
order accurate. This convergence, however, is usually spoiled by the boundary conditions. It is 
also investigated which value of the relaxation parameter minimizes the error. Furthermore, for 
simulations using velocity boundary conditions, an artificial mass increase is often observed. This 
is a consequence of the compressibility of the lattice Boltzmann fluid. We investigate this issue and 
derive an analytic expression for the time-dependence of the fluid density in terms of the Reynolds 
number, Mach number and a geometric factor for the case of a Poiseuille flow through a rectangular 
channel in three dimensions. Comparison of the analytic expression with results of lattice Boltzmann 
simulations shows excellent agreement. 

PACS numbers: 47.ll.-j 



I. INTRODUCTION 

The simulation of fluids using the lattice Boltzmann 
method (LBM) has become a very important research 
field in the last two decades [ij-Q due to its simple and 
parallelizable implementation and various fields of appli- 
cation 

In the literature, the velocity field u of the fluid usu- 
ally is the central observable of interest. More recently, 
however, some researchers are becoming interested in the 
shear stress field a oc Vit, in particular when studying 
blood flow High shear stresses may promote 

blood clotting, an effect with quite diverse consequences 
such as stopping blood flow through injuries (a desired 
effect) or blocking vessels plagued by arteriosclerosis jT6| . 
In the case of diffusive scaling, the time step and the lat- 
tice constant obey At cx Ax^ , and the solution to the lat- 
tice Boltzmann (LB) equation converges to the solution 
to the Navier-Stokes equations with a second order rate, 
i.e., the relative error eu of the velocity scales with Ax^. 
Noting that, within diffusive scaling, the Mach number 
linearly scales with grid resolution. Ma oc Ax, one also 
obtains e„ cx Ma^. 

To our knowledge, a systematic analysis of the conver- 
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gence behavior of the shear stress in the LBM has not 
been conducted yet. It is often stated that the shear 
stress enters the LB equation on the level leading to a 
convergence cx Aa; only. As will be shown in this report, 
this a priori statement underestimates the accuracy of 
the LB simulations. 

Velocity boundary conditions are very useful in the 
LBM. They allow to impose a desired velocity field at the 
borders of the computational domains. This is crucial, 
since the solution to the Navier-Stokes equations crit- 
ically depends on the boundary conditions. One draw- 
back of those boundary conditions is the unphysical mass 
increase, observed under some circumstances. In this pa- 
per, we study this aspect via analytic means. In partic- 
ular, using the compressibility of the LB fiuid, we solve 
the Stokes equation for a rectangular channel subject to 
a Poiseuille fiow and show that the fluid density exponen- 
tially grows with time. The rate of mass increase is found 
to depend on the inverse Reynolds number but scales as 
the third power of the Mach number. We demonstrate 
the validity of analytic results through comparison with 
LB simulations. 

This article consists of two parts. First, we investi- 
gate the convergence of the shear stress in a laminar 3D 
Poiseuille flow. Moreover we address the issue of opti- 
mum parameter choice for the LB simulations. We es- 
pecially search for a reasonable value for the relaxation 
parameter r in view of the shear stress as the relevant 
observable. The second part deals with the global mass 
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increase in the simulations, which arises when velocity 
boundary conditions for both the inlet and outlet of the 
flow are used. This is particularly important for the cor- 
rect computation of the shear stress. 

The paper is organized as follows: In section |TT1 the 
theoretical background is provided. After a short intro- 
duction to the LBM in subsections III Al through FlI CI the 
convergence behavior is discussed in subsection III Dl Ve- 
locity boundary conditions are briefly dealt with in sub- 
section III El and the influence of the initial simulation 
parameters is addressed in subsection III F I The mass in- 
crease is analyzed in section lllll The numerical results 
are presented in section Hvl followed by the conclusions 
in section Ivl 

II. THE LATTICE BOLTZMANN METHOD 
A. The basic concepts 

The LBM has become a very efficient Navier-Stokes 
solver. Its strength is based on its simple coding and its 
locality, which makes it intrinsically parallelizable. Only 
next-neighbor information is required during the evalua- 
tion of the LB equation. 

The LBM can either be regarded as the discretized 
form of the Boltzmann equation or as the extension of 
lattice gas automata [13]. There is one major difference 
between the LBM and other Navier-Stokes solvers: While 
conventional methods directly solve the Navier-Stokes 
equations in terms of the density p and the velocity u 
and approximate the differential equations by finite dif- 
ferences, the LBM introduces a number of q populations 
fi (i = 0, . . . , g — 1) streaming along a regular lattice in 
discrete time steps. Those populations can be regarded 
as mesoscopic particle packets propagating and colliding. 
For the entire LBM, only algebraic manipulations of the 
populations are required. A more detailed presentation 
of the LBM can be found in the literature (e.g., and 
the references therein). 

It is natural to consider all following equations in this 
article as being dimensionless. The physical dimensions 
are recovered by multiplying the lattice quantities by ap- 
propriate combinations of the lattice constant Aa;, the 
time step Ai and the density p given in m, s, and kg/m'^, 
respectively. In lattice units, the time step and lattice 
constant are usually set to unity. 

The evolution of the populations fi is given by the LB 
equation, which takes the form 

f,{x + c,,t+l)- U{x,t) = --{h{x,t) ~ ft'^ix.t)) (1) 

T 

in the Bhatnagar-Gross-Krook (BGK) approximation. 
The equilibrium f^'^ is defined below, r is the relaxation 
time of the LB fluid. At each time step, the populations 
propagate along the q discretized velocity vectors to 
the next neighbors. At those points they collide accord- 
ing to the right-hand-side of Eq. ([T]). In this article, we 



use a 3D model with 15 velocities, designated DSQIS. 
This lattice and the corresponding velocity vectors are 
introduced in The dimensionless relaxation param- 
eter T is connected to the speed of sound Cs and the 
kinematic viscosity v of the fluid by v = Cg(r — 1/2)- In 
lattice units the sound speed takes the value = 1/3. 
In the following we will use Latin indices for the popula- 
tions fi and Greek indices for spatial components, e.g., 
Ua- Throughout the article, the Einstein sum convention 
is used. 

The q populations fi carry the entire information of 
the fluid, and the macroscopic observables like density 
and velocity can directly be recovered by 

q-l q^l 
P=^fl, pu = ^c,fi (2) 

1=0 1=0 

at every fluid lattice node. Also the deviatoric shear 
stress (T can be computed from the populations. We will 
come back to this issue in more detail in subsection III CI 
In the following we drop the summation limits. 

In Eq. ^ the equilibrium populations are given by 

/r = + 3c. -u+^ic,- uf - • • (3) 

This is closely related to the truncated form of the 
Maxwell distribution, which is a very good approxima- 
tion for small Mach numbers. Thus, the LBM is only 
reasonable for small velocities compared to sound speed. 
The q factors Wi are the lattice weights, depending on 
the underlying lattice structure. Their choice ensures 
the isotropy of the fluid, a necessity to solve the Navier- 
Stokes equations asymptotically. For the -DSQIS lattice, 
they are also given in [3]. 

Beside the constraint of small velocities, the LBM suf- 
fers from another shortcoming: the incompressibility of 
the fluid gives way to a quasi incompressibility indicated 
by the ideal gas equation of state, p — c^p. As will be 
exposed in section IIIII this is the reason for the mass 
increase when using velocity boundary conditions. 

In order to show that the LB equation is equivalent to 
the incompressible Navier-Stokes equations in the limit 
Ax — > and Ma 0, one performs a Chapman-Enskog 
analysis, which will be briefly presented in the following 
subsection. 



B. Chapman-Enskog analysis 

The idea behind the Chapman-Enskog analysis applied 
to the LB equation is that different physical phenomena 
happen on different time scales. While the advection 
of the fluid is the fastest process, the diffusion of mass, 
momentum and energy happens on a slower time scale. 
For this reason, the time derivative is usually split into 
two parts. One can also extend this approach and take 
into account more time scales, 

dt = edto+e^dt,+e^dt,+0(e^). (4) 
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Since the spatial variations of all processes are of the 
same order, the gradient is not decomposed, V — eVi. 
The populations are expanded about the equilibrium, 



(5) 



where fj^^^ — f°^. The parameter e <C 1 can be identified 
as the Knudsen number (the ratio of mean free path to 
a characteristic length, usually the size of an obstacle or 
the entire system). The Taylor expanded LB equation 
[Eq. (P)] takes the form 

°° 1 1 

^ - (ft + c, . V)" Mx, t) = - {Mx, t) - ft\x, t)) . 

Tl. T 

n—1 

(6) 



The expansions for ft, V and fi are inserted into Eq. 

and the terms are sorted by the powers of e. The 
e-equation reads 



(ft„ + c, . V)/, 



(0) 



the e^-equation 



(1) 



+ c.-V)/« 

c..v)Vr 



f(2) 



(7) 



(8) 



(2) 



(0) 



■to 



c.-VffP 



f(3) 



(9) 



For recovering the Navier-Stokes equations in terms of 
p and u, the first and second order equations ([7]) and 
^ and the same equations multiplied by Cj are summed 
over i. One can show that this approach yields the four 
macroscopic equations 

dt,p + v ■(pu) = o, ft„(pM) + v-n(") = 0, (10) 



ft,p = o, ft,(ptt)+ (^i-^j v-n(i) = 0. (11) 

The tensor 11 has the components 

Ha^ = ^ CiaCtpfi , 
i 

i 



The incompressible Navier-Stokes equations read 

V-M = 0, pdtu + pV -{uu) = -Vp+2i^pV ■ S . (15) 

In the incompressible limit (Ma — > and Ax -> 0), Eqs. 
(1131) and (jlSp are equivalent, if we assume as equation of 
state p = c^p for the LB fluid. This means that the LBM 
asymptotically solves the incompressible Navier-Stokes 
equations. For a detailed Chapman-Enskog analysis we 
refer to [3]. It is also possible to compute the macro- 
scopic equations of higher orders in e. In 19] the result 
for is presented. 

As a further comment, we point to non- isothermal lat- 
tice models. The standard LBM, as presented in this 
article, does not satisfy the Galilean invariance princi- 
ple when considering moments of order higher than two. 
For the numerically stable modeling of thermal systems, 
which depend on higher moments, extended lattice struc- 
tures have to be considered I20l-[23l. 



Cicc denotes the a-component of the velocity vector c^. 
IIq^ will be related to the macroscopic momentum flux 
tensor in subsection III CI The combined form of the two 
time scales to and ti finally reads 



ftp + V • (pu) = , 
ft(pit) + V • (puu) = -c^Vp + 2i^V • (pS) . 



(13) 



ly — c'^{t — 1/2) is the kinematic viscosity of the LB fluid. 
The shear rate tensor S has the components 
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(14) 



C. Shear stress in the lattice Boltzmann method 

Fluids (liquids and gases alike) are described by the in- 
compressible Navier-Stokes equations, cf. Eq. p^ . in the 
limit of small Knudsen and Mach numbers. Those equa- 
tions can also be written in the form pftwa = —dpMap 
with the total momentum flux 



Maf} ^ pSaf3 + pUaUp - Gaf} 



(16) 



where 5ai3 is the usual Kronecker symbol. The first term 
on the right-hand-side of Eq. (fT6|) is the isotropic pres- 
sure and the second one the momentum transfer due to 
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mass transport. The third term describes the momen- 
tum diffusion due to the viscosity of the fluid, ft is called 
the deviatoric shear stress. It only vanishes if the veloc- 
ity gradients are all zero or if the entire fluid is rotating 
with a spatially constant frequency (as to ensure that 
no relative motion occurs within the fluid). In a vis- 
cous fluid, momentum diffuses from regions with large to 
regions with small momentum. If the fluid is incompress- 
ible, this deviatoric shear stress tensor, from now on only 
called shear stress tensor, takes a simple form. 



Tap = V [daUp + dpUa) 



(17) 



rj = pv is the dynamic viscosity. For Newtonian fluids, 
the shear rate tensor S Eq. (|14p and the shear stress 
tensor cr Eq. (|17p are related through 2pvS ~ cr. Both 
tensors are symmetric and traceless. The latter prop- 
erty holds as long as the fluid is incompressible, since 
Tr S' = V • It. It is known that some fluids show a shear- 
thinning or shear-thickening behavior, i.e., the viscosity 
is a function of the shear rate (e.g., blood or honey). If 
the fluid is Newtonian, viscosity is not a function of the 
shear rate, and 5' and cr differ only by a constant fac- 
tor. Since we restrict ourselves to the latter case, unless 
otherwise stated, the terms 'shear rate' and 'shear stress' 
can be used synonymously. 

It follows from the Chapman-Enskog expansion that 
the momentum flux tensor M from the Navier-Stokes 
equations can be approximated by the second moments 
of the LB populations. The equilibrium parts of the pop- 
ulations lead to the pressure and the convective term in 
the momentum flux tensor, 

n^°j = ^ CiaCipff'^ = C^spSaf} + pUaUjj . (18) 
i 

The shear stress is provided by 

i 

= [da{pUfj) + dp{pUa) + dtoipUaUp)] (19) 

r 

Obviously, there appear some correction terms in Eq. 

which need to be discussed. We will come back 
to this issue in subsection IIIDI The total Navier-Stokes 
momentum flux can be approximated by {p = c^p) 



Map = p5aP+pUaUp-2pvSaP ~ 11 



(0) 

ap 



1 _ _ 1 n^i) 



(20) 

where ly = c'^{t — 1/2) has been used. 

At this point a difficulty arises. It is easy to sepa- 
rate the equilibrium populations f^^"^ = f^'^ from the 
remaining corrections, the non-equilibrium populations 

jnoq _ _ joq pOSSlblc tO COmputC 

the first order populations f^^^ as stand-alone quanti- 
ties. This means that /"'''^ has to be used instead of f^^K 



Obviously, this additional approximation only concerns 
the shear stress tensor. The remaining macroscopic vari- 
ables (density, velocity) are not inffuenced. Sometimes 
the tracelessness of the shear stress tensor is enforced 
explicitly by writing 



(Tap 



1 

27 



E 



^ia^-i^ ^ ' 



(21) 

where D is the number of spatial dimensions. Due to the 
approximation /^^^^ — >■ and especially the influence 

(21 

of /j , it is generally expected that the LBM result for 
the shear stress converges with a first order rate in the 
case of diffusive scaling {At cx Acc^) — in contrast to 
the velocity, which has a second order accuracy. This 
convergence behavior is further discussed in subsection 
IIIDI and tested in section IIVI 

Up to this point we have not mentioned the special role 
of the shear stress in the LB simulations: Spatial deriva- 
tives on a lattice are usually computed by using finite 
difference (FD) schemes. The advantage of the LBM is 
that the shear stress — although related to the gradient 
of the fluid velocity — can directly be assessed locally, 
i.e., on each individual lattice node, and independently 
of the velocity. No information of neighbors is required, 
and the implementation is straightforward. This is the 
strength of the LBM. Nevertheless, in order to have a 
comparison, we have examined in section IIVI both the 
shear stress obtained from the local approach given in 
Eq. (PT|) and a non-local finite difference method. All fi- 
nite differences in this paper are of second-order accuracy 
and central, e.g., 

du^ f ^ Ua:{x,y + l,z) - u^{x,y - l,z) 

= . (22) 



D. Convergence of the shear stress 

Knowing the Navier-Stokes equations and the macro- 
scopic limit of the LB equation, the basic idea is to es- 
timate the error of the solution obtained by the LBM, 
as compared to the Navier-Stokes solution. As also men- 
tioned in the introduction, in the diffusive scaling, the 
relative error of the velocity is proportional to Ma^ (see 
e.g., [131 and references therein). The same is true for 
the density fluctuations, i.e., the compressibility of the 
LB fluid also scales with Ma^ . The convergence behavior 
of the shear stress, on the other hand, is not obvious and 
has to be deduced from the LB equation. 

Since the equilibrium populations can be computed ex- 
actly (in terms of the Maxwell truncated approximation), 
the equilibrium part 11^°^ Eq. contains only com- 
pressibility and truncation errors of the same order as 
the velocity and the density. Higher order terms of fi do 
not play a role, since they do not appear in the equation. 
In the diffusive scaling, the relative velocity and density 
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errors are c>c Ma^ and Cp cx Ma^ and, for this reason, 
the relative error of n|^'^^ also scales with Ma^ . 

Things are different for the shear stress Il^^j in Eq. 
(Uni). We write the LBM shear stress as the Navier-Stokes 
shear stress plus an error term, 

nS ^-\p i^cUf! + dpuo,) + Eni . (23) 



It can be shown [18| that this error term reads 

(24) 

The point here is that not only the errors of p and Ua 
account for the total deviation of the shear stress, but 
also the additional terms in Eni . For the relative error 
it follows that eni = ii'ni /nj^^'j cx Ma^, since the three 

relations Uq oc Ma, da oc Ma, and daP oc Ma"^ hold. 
The first statement is obvious. The second one can be 
understood by recalling the diffusive scaling Aa; oc Ma. 
When the grid is refined and Aa; decreased, all spatial 
derivatives in lattice units are also decreased by the same 
rate, i.e., da oc 1/N oc Ax. The last statement follows 
from the second one and the fact that the LB fluid com- 
pressibility scales with Ma^, i.e., the fluctuations of the 
density p with respect to the mean density p behave like 
{p- p)/p(x Ma^. 

This indicates that the shear stress would converge 
with a second order rate, if the populations f^^^ were 
known and used for its computation. As we have stated 
before, it is not possible to evaluate f^^\ In principle, the 

populations f^^^ can be computed from the equilibrium 
via Eq. ([7]). This is not practical, because information 
from next neighbors and different time steps would be 
necessary in order to evaluate this equation. The LBM 
obviously would loose the advantage of locality, and the 
shear stress could directly be computed using a finite dif- 
ference scheme for the equation aai3 — vp{daUf} + dpUa) 
in the first place, without bothering computing f^^^ at 
first and 11^^"'^ afterwards. Hence, the first order terms 
f^^^ are substituted by the non-equilibrium populations 
in the simulations and is computed instead of n^"^j . 

As a consequence, the influence of the higher order cor- 
rections /j^^^'' plays a role. Fortunately, the higher order 
populations become less important with increasing order 

(2) 

of e. Hence, we focus only on the influence of Jl and 
check whether its contribution may destroy the second 
order convergence of the shear stress. One can write 



and thus 



K7 - - gP {daU0 + dpUa) + Eui + E, 



(25) 



E^ is the error due to the higher order corrections, espe- 

(2) 

cially /j . Thus, we are interested in the convergence 



behavior of H^^^^ = J^i'^ia'^i/sfi^' ■ From the Chapman 
Enskog analysis we know the population fj^"^^ , cf . Eq. ([5]) 



(1) 



13 ^ "l-^-afSj 



(26) 

where Rais-y — X^i CiaCipCi^fi is the third moment. From 
Eq. ^ it follows 



(27) 



with the fourth moment P, 



(0) 



(0) 



Since the standard equilibrium Eq. ([3]) is truncated af- 
ter the second order in the velocity, all velocity moments 
larger than the second order are not physically well de- 
fined. Indeed, in the framework of the truncated equilib- 
rium, all equilibrium moments larger than two scale with 
pu^ oc Ma^. 

At this point we have to admit that it is a tedious 
task to exactly compute 11^^^ . Actually, this is not neces- 
sary, because the focus is put on the convergence behavior 
and thus the scaling with the Mach number. The follow- 
ing considerations are straightforward if one recalls that 
Ua oc Ma, da oc Ma, and daP oc Ma'^: From the above 
observations we conclude that all terms in Eq. (|26p are of 
order Ma''. Additionally, Eqs. (fTTjl . ([1^ . and ^ 

have been used. This means that both error terms Eni 
and _Ee are of fourth order in the Mach number, whereas 
n^^^ itself is second order. As a consequence, the rela- 
tive error of the deviatoric shear stress, as computed from 
Eq. , is quadratic — similar to the errors of the veloc- 
ity and the density. This is a quite encouraging result. 
However, as will be shown below, the second order con- 
vergence is very sensitive to additional error sources. In 
particular, our simulations clearly show that the conver- 
gence behavior strongly depends on the type of boundary 
conditions. 

In subsection llV Al our numerical results are presented 
and compared to the above qualitative findings. 



E. Geometry and boundary conditions 

For the simulations in this article we have chosen a 
simple 3D geometry, the laminar 3D Poiscuillc flow in a 
channel with rectangular cross section. The fluid enters 
the numerical grid at a; = 1 and exits at a: = N^- The 
channel has height H (along z-axis), width W (along y- 
axis), and length L = [N^ — l)Aa;. The Reynolds and 
Mach numbers are 



Hu ^ ^ u 
Re = , Ma = — 



(28) 



u is the center velocity of the flow. 

The velocity field for an infinitely long channel with 
rectangular cross section in 3D is known analytically. 
The Navier-Stokes equations reduce to the Poisson-like 
equation •q/S.Ux{yTz) = dp/dx = const. It can be solved 
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by expanding the solution in basis functions. One possi- 
ble solution is 



Ux{y-,Z) = — 



z{H ~ z) 



^ ^ cosh(n7r(y - W/2)/H) ^ (7mz\ 



■E 



n3 cosh(n7rT/F/2i?) 



H I 



(29) 



with the origin of the coordinate system at one corner of 
the channel. The pressure gradient has been eliminated 
in favor of the center velocity Ux, dp/dx = —2r]Ux/^, 
where E is a normalization factor depending on the ge- 
ometry. For the case of the rectangular channel described 
above, 



^2 8H^^ 1 



sin(n7r/2) 



4 TT^ ^ cosh(n7rH//2i7) 
Another equivalent solution 25] is 



(30) 



^ or, odd -I 



sin sin 2f2 



n,m ^ ' ' \ I J 



(31) 



with the same value of E. 

Both solutions have advantages and disadvantages. 
The first solution has only one expansion coefficient and 
converges faster. However, while convergence is fast in 
the region near y — W/2, the no-slip condition at y = 
and y — W is poorly reproduced. The second solution, on 
the other hand, shows very good convergence near the en- 
tire border, since the basis functions vanish there. How- 
ever, two expansion coefficients are necessary. Moreover, 
the convergence is slower when approaching the center 
of the channel. Therefore, when computing the analytic 
solution, we divide the system into an inner part, where 
Eq. (jH E[) is used, and the wall region, where the exact 
solution is computed via Eq. (|3ip . The shear rate tensor 
S can directly be calculated from both solutions. 

The accuracy of the simulation is quantified using the 
common L2-norm for the relative error 



(32) 



where the sum takes into account all lattice nodes in 
the numerical grid. The error of the two relevant shear 
components Sxy and Sxz are defined similarly. 

The no-slip condition at the side walls is realized by a 
fuUway bounce-back rule. The walls are located halfway 
between the fluid and the wall nodes, and the accuracy is 
of second order [2^ . For the inlet and outlet we have used 
periodic boundary conditions in one part of the simula- 
tions and velocity boundary conditions in the other part. 
In the case of periodic boundary conditions, the flow is 
driven by a constant and homogeneous body force [27| . 



This force can be implemented by adding Wi Ci ■ F/cj. to 
the right-hand-side of Eq. ([T]), where F is the body force 
pointing in x-direction. 

We have examined both the velocity boundary condi- 
tions proposed by Skordos ^] and by Latt ^29J. The 
velocity profiles on the inlet and outlet are fully devel- 
oped. Implementing velocity boundary conditions is not 
as straightforward as using the bounce-back rule. The 
underlying idea is the following: If the fluid velocity u 
is given on a straight wall, the populations fi have to be 
computed from u. This is in general a problem, since the 
number of equations available is smaller than the number 
of unknown populations. For this reason, some closure 
relations are assumed, which complete the set of equa- 
tions. There are various possibilities to do this. A good 
review of the velocity boundary conditions can be found 
in [29|. Using an explicit method, it boils down to the 
equation 



(1) 



(1) _ W^Tp 



(33) 



with Qia0 — CiaCii3 — c^Sa/s- FoUowing this general ap- 
proach, the boundary conditions also show a second or- 
der convergence. This is important in order not to spoil 
the quality of the overall simulation. It depends on the 

method how f^^'' is estimated. Skordos uses a finite dif- 
ference method for the velocity gradient, while Latt ap- 
proximates the unknown values of fj^'^'^ with the help of 
the known ones. The density p at the inlet and outlet 
can be recovered by 

p=-^{2p+ + po), (34) 

where uj_ is the projection of the boundary velocity u 
on the boundary normal, pointing outside the numeri- 
cal grid, i.e., u±_ is negative at the inlet and positive at 
the outlet. This way, the correct pressure gradient is 
recovered. p+ is the sum of all populations leaving the 
numerical grid (those populations are known), and po 
is the sum of all populations streaming parallel to the 
boundary plane (those are also known). 



F. Influence of the relaxation parameter 

The accuracy of a LB simulation depends on the LBM 
itself, the initial conditions, the boundary conditions and 
also on the simulation parameters. There are four rele- 
vant quantities which have to be set up for every simu- 
lation: The Reynolds number Re which is fixed by the 
physical system, the Mach number Ma of the simulation, 
the dimensionless relaxation parameter r and the lattice 
constant Ax or, equivalently, the number N of fluid lat- 
tice nodes along one axis (we use N^). Using Reynolds 
and Mach numbers from Eq. (pS)) yields an important 
relation between those four dimensionless parameters. 



Ma 
Re 



1 r- 



(35) 
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Assuming that the Reynolds number is fixed, there are 
two degrees of freedom in the choice of the parameters. 
Although the Mach number in reality is uniquely defined 
by the speed of sound, it can be varied in the LBM. As 
long as the Mach number is small, Ma < 0(0.1), it does 
not influence the physics of the system. This way, in a 
simulation the value of Ma can be adjusted in order to 
increase the accuracy or shorten the computing time. 

The second order convergence of the LB equation can 
only be realized, if the Mach number and the lattice con- 
stant are decreased simultaneously. The diffusive scaling 
fulfills this requirement [13, 113], and it is ensured that the 
incompressible Navier-Stokes equations are recovered in 
the limit of Ax 0. The disadvantage of this method is 
that by reducing the Mach number the required number 
of time steps is increased (recall that At (x Ax^). As a 
consequence, halving the lattice constant and the Mach 
number simultaneously reduces the error by a factor of 2^ 
but increases the simulation time by a factor of 2^ (2D) 
or 2^ (3D). 

Equation (|35p plays an important role by the selection 
of simulation parameters. In order to see this, we note 
that the choice of the Reynolds number and the Mach 
number is closely linked to the physical problem of inter- 
est. In the diffusive scaling, on the other hand, the Mach 
number uniquely sets the lattice resolution Ace and the 
time step At (recall that Ace cx Ma and At oc Ma^). 
Thus, for fixed values of Re and Ma, the only freedom is 
in a variation of the relaxation parameter r or, equiva- 
lently, the number of grid points Nz. 

In many cases of interest, however, the Reynolds num- 
ber and the Mach number can be varied in a reasonable 
range without changing the basic physics of the problem. 
When studying problems, where sound propagation does 
not play a role, it is often possible to increase the Mach 
number as long as density fluctuations remain negligibly 
small. This brings the advantage of using larger grid reso- 
lutions and hence significantly reducing the computation 
time. A similar philosophy also applies to the choice of 
the Reynolds number: When studying laminar flow prob- 
lems. Re can be increased as long as velocity fluctuations 
remain negligibly small. A more complete description of 
these and similar ideas can be found in [3]| . 

Therefore, the question arises which initial choice of 
the four parameters appearing in Eq.[3niis optimal, mean- 
ing which choice leads to the smallest possible error at 
constant simulation time. There is a very detailed anal- 
ysis of the LB error in 24]. The authors point out that 
the accuracy of the LBM strongly depends on the choice 
of the relaxation parameter r. This can also be seen in 
[l9j . The important point is that the error has a min- 
imum for T in a region between 0.8 and 1.0, depending 
on the Reynolds number. Setting t too close to 1/2 or 
larger than 1 results in avoidable inaccuracies of the sim- 
ulations. 

In order to verify the results presented in [13] and to 
see whether the velocity error and the shear stress error 
behave similarly we have set up a series of simulation 



according to the following paradigm: 

1. Set Re to the desired value. 

2. Choose a reasonable value for Ma. 

3. There is only one free parameter left. Set up a 
series of simulations for different grid sizes N. This 
automatically fixes r according to Eq. ([55]) . It is 
wise to start for t around 0.9. 

4. Calculate the relative errors and es of the sim- 
ulations. The curves eu{T) and es(T) both feature 
a minimum, which can be found by varying N and 
thus T. 

5. The minimum value Topt is the best choice for the 
simulations. It should be kept fixed, if the resolu- 
tion is to be changed. This directly leads to the 
diffusive scaling and a second order convergence. 

There are some remarks: In the simulations, the min- 
imum value of T is not necessarily the minimum value of 
the LBM, because the accuracy of the boundary condi- 
tions also depends on r and — in general — in another 
way than the LBM. The bounce-back rule for example 
introduces a slip velocity at the walls, which depends on 
the relaxation parameter Additionally, the qual- 

ity of the velocity boundary conditions depends on the 
Mach number. For this reason, the location of the opti- 
mum value Topt is also a function of Ma for larger Mach 
numbers. 

We present our results for Topt in subsection IIVBI 

III. MASS INCREASE 

Some wall boundary conditions do not obey local mass 
conservation. This shortcoming can be removed by ap- 
plying modified boundary conditions (s^ . [ssj . But even 
if the mass is conserved at the walls locally, it is pos- 
sible that the total mass of the numerical grid changes 
with time . This effect is caused by inlet and out- 

let velocity boundary conditions. Due to the equation 
of state of the lattice fluid, p — clp, a pressure gradient 
is equivalent to a density gradient, leading to a larger 
inlet than outlet flux, and mass is accumulating in the 
system. There are possibilities to oppose this effect, e.g., 
the use of velocity inlet / pressure outlet boundary condi- 
tions or incompressible LB models, both for steady and 
unsteady flows [l^, H^. Mass increase is fully absent 
also in the case of a body force driven flow with periodic 
boundary conditions. However, we are interested in the 
consequences of the mass increase, if it is not avoided in 
the first place. 

The reader should be reminded that L, W and H are 
the length, width and height of the channel, whereas N^, 
Ny and Nz are the number of fluid lattice nodes along the 
length, width and height of the channel, respectively. In 
lattice units, obviously L = N^ — l, W = Ny and H = Nz 
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hold, as long as inlet / outlet are defined by velocity 
boundary conditions and the side walls by the bounce- 
back rule. The side walls themselves are not contained 
in Ny and A^^, so the total lattice — including fluid and 
walls — has dimensions x [Ny + 2) x {Nz + 2). 

For a 3D rectangular Poiseuille flow, the velocity pro- 
files at the inlet and the outlet have to be equal, u — 
Mill = '"out; resulting in a net total mass increase per 
time step 



M 



AM 



-ApAu > 0, 



(36) 



where Ap is the density difference between the outlet and 
the inlet (and thus negative), A is the cross section area 
and u the mean velocity on the cross sections. Note that 
u and u are proportional to the Mach number and that 
tt/u is constant for a given channel geometry. In a 2D 
Poiseuille flow one flnds u/u — 3/2. In three dimensions, 
on the other hand, this ratio depends on the aspect ratio 
W/H of the cross section. 

In an ideal LB simulation of Poiseuille flow along the x- 
axis, the density is only a function of x. Furthermore only 
the velocity component Ux does not vanish. It is constant 
along the x-axis and depends on y and z. In this case, 
the Stokes equation for the LB fluid reads (assuming that 
steady state is reached, i.e., dtu ~ 0) 



dlnp 
dx 



const , 



(37) 



whose integrated form reads 



Qvu 

Pa cxp [ — —X 



(38) 



E is the normalization factor from Eq. (150]) . The pressure 
difference between inlet {x — 1) and outlet {x = Nx) is 



Ap = po 



cxp ( — ^iVa; I - cxp I — — 



(39) 



The inlet and outlet densities are given only implicitly by 
Eq. It is, therefore, not possible to set the velocity 

profiles and the densities on the inlet / outlet simultane- 
ously. For small values of GvuNx/T,, the density gradient 
can be approximated by the expression from subsection 
UlEl -QpovuL/Y., i.e.. 



Ap 

T 



(40) 



This is valid whenever the Mach number is small and the 
density near pg everywhere. 

The total mass in the numerical grid can be estimated 

by 

Af w / dV p = A / dxp 

Jl/2 



-Apo—exp — — 
bvu \ 1^ J 



exp 



-—Nx 



Putting everything together results in the expression 



M 
M 



Qvuu 1 — exp 



K S 



{N. - 1)) 



exp 



{^Nx 



■ exp 



&VUU 



E 

(42) 

The approximation in the second step is justified if Nx 
and E are sufficiently large, which is always the case for 
an appropriate lattice resolution. 

Since the geometry and the Reynolds number shall be 
fixed here, the Mach number can only be adjusted by 
changing the kinematic viscosity u according to oc Ma 
so that finally 



M{Nt) = Mo exp {CgMs?Nt/Re) 



(43) 



The mass after Nt time steps depends on the Mach and 
Reynolds numbers and also on a geometry factor 



2 Hu 



1 



x/3 S 



H 



(44) 



which is uniquely fixed by the shape of the channel. The 
Mach number enters the equation with a power of 3. It 
is an observation that usually CgMa'^/Re <C 1, but for 
a large number Nt of time steps the mass increase may 
enter its exponential regime. When examining the valid- 
ity of Eq. (|43p . one must pay attention to the following 
aspects: 



1, 



The equation is not valid at the very beginning of 
the simulations, but only when a 'steady state' has 
been reached in which the pressure gradient is fully 
developed. Here, steady state means that the ve- 
locity does not change in time, although the density 
is not constant. For this reason, the first iteration 
steps are neglected when comparing simulation re- 
sults with the analytic expression Eq. 

2. Modifying the Mach number over a wide range, but 
keeping Re and N fixed requires changing t over a 
wide range. One has to be careful that t stays in 
a region, where the boundary conditions and the 
LBM work reliably. 

We will show in section lTVl that once the above considera- 
tions are taken into account, Eq. produces excellent 
results. 

Another question is how a change of the simulation pa- 
rameters influences the mass increase. If Re is changed, 
this has to be compensated by varying the Mach num- 
ber, the grid resolution or the kinematic viscosity or 
any combination of those. Suppose the following in- 
dependent transformations: Re a Re, Ma — ?> 6 Ma 
and cN^. Since Cg cx 1/7V^, Nt oc N^/iy, 

Xjv — Re/{Nzu) and 1/u oc 1/Ma, one flnds that 
CgMa^TVt h'^CgMs?Nt. For fixed aspect ratio of the 
channel dimensions, this means that the mass increase is 
directly related to the Mach number and is independent 
of the other parameters. The compressibility (controlled 
by the Mach number) and the mass increase are tightly 
(41) connected. 
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The mass increase actually raises a problem: The fluid 
velocity within the LBM is defined in such a way that a 
global variation of the mass does not change its value, as 
long as the pressure gradient is unchanged and the distri- 
bution of the populations is affected only by a velocity- 
independent prefactor, cf. Eq. The definition of the 
shear stress in Eq. (|2ip on the other hand is sensitive to 
the global mass, since cr oc p. If the non-physical mass 
increase is significant, the shear stress tensor must be 
corrected for this effect. 

There are in principal two options for this correction. 
First, one may take measures to compensate the mass 
increase in the first place. Second, one may allow the 
mass increase, but make sure that all observables remain 
unaffected. We assume that the shear rate tensor Eq. 

is not affected by the mass increase, since it does 
not depend on the density. Hence it is reasonable to 
calculate this quantity instead of the shear stress tensor. 
Due to 2vpS = cr, one finds, starting from Eq. ([21]), 



Ma 



0.16 



0.08 



0.04 



0.02 



flicq 



(45) 



In order to check whether this approach is reasonable, 
we have computed the relative error es of the shear rate 
as a function of time in a 'steady state' simulation with 
significant mass increase. The results are discussed in 
subsection IIVCI 



IV. SIMULATIONS AND RESULTS 

We investigate the velocity and shear stress error be- 
havior in a laminar 3D Poiseuille fiow using both periodic 
boundary conditions with body force and velocity bound- 
ary conditions. We are interested in the convergence of 
the numerical results towards the analytical solutions. 
Moreover, we determine the optimum choice of the re- 
laxation parameter, depending on whether the velocity 
or the shear stress is the relevant observable. Those re- 
sults are presented in subsections IIV Al and IIV Bl 

The computing time of a LB simulation can be re- 
duced if the Mach number is increased. The reason 
is that less time steps are necessary to reach the final 
state of the flow. This is very important whenever large 
Reynolds numbers are given, because simulations with 
high Reynolds number require a large numerical grid. 
However, the increase of the Mach number results in 
a larger numerical error due to compressibility effects 
and, in combination with velocity boundary conditions, 
it leads to the violation of mass conservation in the sys- 
tem. This mass increase in a steady Poiseuille flow is 
examined qualitatively and quantitatively in subsection 
HVCl 

The presented benchmark analyses are not intended to 
and cannot cover the above mentioned aspects in every 
detail. Rather, our aim is a fundamental understanding 
of the discussed observations. Nevertheless, it would be 
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Figure 1. Relative L2 errors of the velocity Ux and the shear 
stress component Sxy (both local and FD) as functions of the 
dimensionless system size A'': The velocity slope is close to —2 
for all boundary conditions, whereas the shear stress slope is 
between —1.4 and —1.9 [local via Eq. (|45|) ] and around —1.5 
(FD). The slope values are collected in table H] 



an interesting task for future work to extend the above 
studies to other situations such as complex geometries 
as well as non-steady flows, issues not considered in this 
work. 
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Figure 2. Behavior of the relative shear error es(r) [left column: from Eq. (|45p . right column: from FD] for Reynolds numbers 
1 (top), 10 (middle), 100 (bottom). 



A. Convergence 

In subsection IIIDl we have argued that the conver- 
gence of the shear stress should in principle have a Ax~^- 
or Ma~^-behavior. This prediction has been tested with 
a series of simulations at fixed Re = 12.5 and r « 0.855. 



The size of the numerical boxes has been 16 x 18 x 18 up 
to 128 X 130 X 130. The errors of the velocity and the 
shear component Sxy are shown in figure [T] as a function 
of the lattice resolution. 

There are three relevant observations. First, the rela- 
tive error e,t of the velocity is significantly smaller than 
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that of the shear rate. This is no surprise, since there are 
many error terms in Eq. (|25p . 

The second and most interesting observation is the 
slope of the convergence. For the velocity, it is found to 
be very close to —2 for all boundary conditions, see table 
m It shows that the LBM and the boundary conditions 
are second order accurate with respect to the velocity. 



Table I. Slope values of the velocity Ux and shear rate Sxy 
errors for the three boundary conditions (BC). The corre- 
sponding figures are presented in Fig. [T] 



BC 


Ux 


Sxy (local) 


Sxy (FD) 


periodic 


-1.99 


-1.88 


-1.49 


Latt 


-1.94 


-1.40 


-1.52 


Skordos 


"1.97 


-1.62 


-1.56 



This has been shown many times in the literature. For 
the shear rate Eq. (|45p on the other hand, we have found 
slopes strongly depending on the boundary conditions. 
The slope for the periodic boundary condition is —1.88 
and nearly reaches the predicted value of —2. The ve- 
locity boundary conditions have flatter slopes, cf. table 
m but are clearly better than first order convergence. 

Our interpretation is that indeed the shear stress is a 
second-order accurate observable, but the convergence is 
detrimentally influenced if velocity boundary conditions 
are used. Since the implementations of those conditions 
usually are only self-consistent on the e- and e^-level, ad- 
ditional errors appear affecting the shear stress. 

Another interesting point is that the error of the shear 
stress obtained from the finite difference scheme is not 
significantly smaller than that of the local scheme. The 
slopes for the finite difference shear rate are around —1.5 
for all investigated boundary conditions. 

Let us briefly compare the simulations: The velocity 
boundary conditions yield a much better accuracy for 
the velocity field than the periodic boundary condition 
(nearly one order of magnitude). For the shear stress, 
things are different: The periodic boundary condition 
with a body force gives rise to the highest accuracy and 
the best convergence, but Skordos' method is not signif- 
icantly worse. Latt's velocity boundary condition is less 
accurate, but a 1% accuracy can be found even at small 
resolutions. Additional information on the accuracy of 
the shear stress is provided in Sec. IIVBI 

Some remarks have to be made here. Steady Poiseuille 
flow has the simplest possible geometry. The error of the 
shear stress is given by the terms Eui and i?c in Eq. (j25p . 
It is expected that those quantities can be much larger in 
arbitrary flow geometries and non-steady flows, although 
the second order convergence should be preserved. The 
small relative errors es — 0.0001 . . . 0.01 which have been 
recovered in these benchmark tests may be a consequence 
of the high symmetry of the Poiseuille flow. For this 
reason, more benchmark tests for other, more complex 
geometries, for which the analytical solution is known, 
should be performed. Practically, the second order rate is 
not assumed to be recovered as long as velocity boundary 
conditions are used. 

It is an open question whether there exists a global slip 
velocity in our simulations, which could be subtracted 
from the velocity profiles. If the slip velocity is constant 
over the cross section of the flow, i.e., if the slip velocity is 
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the same at edges and corners, it would be important to 
identify its value and subtract it from the velocity field. 
Such a correction is not included in the present work. 



B. Initial parameters 

For the Reynolds numbers 1, 10 and 100, the optimum 
values of the relaxation parameter r have been investi- 
gated according to the procedure described in subsection 
IIIFI The error curves e(T) are depicted in figures [2] and [3l 
and the optimum values along with the minimum errors 
are given in table HIl 

From a theoretical point of view, Topt depends on the 
value of the Reynolds number [23] ■ The analysis in (2^ 1 
however does only account for the error due to the LBM. 
The behavior of the boundary conditions is not included. 
Those effects are clearly visible in our results. 

For the velocity in figure [3] we have found that veloc- 
ity boundary conditions are in principle more accurate 
than a body force driven flow. At larger Reynolds num- 
bers, however, when the Mach number is increased in 
order to shorten the simulation time, the velocity bound- 
ary conditions become more inaccurate. But even at the 
comparably large value Ma — 0.20, the accuracy of all 
three methods is similar. It strikes the reader that the 
relaxation parameter should be somewhere between 0.9 
and 1.0 in order to minimize the error Cu- The exact 
value depends on the boundary condition and it will most 
likely also depend on the geometry and the Mach number 
(which we have not tested) . In table |TT1 the locations of 
the minima are given. 

The behavior of the relative shear error es is presented 
in figure [2] Both methods (shear locally from Eq. ([45]) 
and a finite difference scheme) can directly be compared. 
The local shear error is smallest for the body force driven 
flow, especially at larger Mach numbers. This observa- 
tion matches the previous discussion in subsection IIV Al 
where it is argued that the velocity boundary conditions 
give rise to additional errors. Both velocity boundary 
conditions behave differently. Skordos' method is more 
accurate than Latt's scheme. While the former boundary 
condition produces errors similar to the body force driven 
flow at Ma — 0.04, the latter method begins to become 
inaccurate. Finally, at Ma — 0.20, the error obtained 
from Latt's method is nearly two orders of magnitude 
larger than that of the body force driven flow, and also 
Skordos' methods becomes less accurate. The position of 
the minimum error is at r < 0.9 for the reasonable simu- 
lations and becomes larger when the boundary conditions 
are inaccurate. The reason is that large r is equivalent 
to high grid resolution (or small Ax) thus reducing the 
error related to a higher Mach number. From the plots 
in figure [2] we learn that the accuracy of the shear suf- 
fers from velocity boundary conditions at relatively large 
Mach numbers. 

This shortcoming can be avoided if a finite difference 
scheme for the shear is used. Here, the relative error is 



tightly connected to the error of the velocity, since the 
shear is interpolated from the velocity information. It 
may be surprising that the minimum error of the shear 
from a finite difference scheme is not significantly smaller 
than that of the shear obtained from Eq. (|45p . So there is 
no a priori advantage of using a finite difference scheme, 
but — at larger Mach numbers — the finite difference 
scheme is clearly superior, at least for the velocity bound- 
ary conditions. 

In conclusion we briefly summarize the most important 
observations: 

1. The minimum error is always located somewhere 
around the interval [0.8, 1.1] for the relaxation pa- 
rameter r, depending on which quantity is com- 
puted and which method is used for this. It is 
therefore very reasonable to use a relaxation pa- 
rameter out of this interval for simulations, r ~ 0.9 
is a compromise for most simulations. 

2. For the velocity field, velocity boundary conditions 
are more accurate up to Ma = 0.20. 

3. One has to act with caution, if the error of the 
shear shall be minimized. Velocity boundary condi- 
tions become increasingly inaccurate for large Mach 
numbers, if Eq. (|45p is used. A finite difference 
scheme overcomes this disadvantage, but the im- 
plementation is more difficult. 

4. Keeping the Mach number small in general is 
recommended, but not always practical in large 
Reynolds number simulations. If a simulation re- 
quires a large Reynolds number, one should de- 
crease T and increase Ma reasonably, keeping in 
mind which observable {u or cr) is of most interest. 



C. Mass increase 

A series of simulations with Reynolds numbers 5, 10, 
15 and 20 and Mach numbers between 0.05 and 0.3 have 
been performed in a channel with dimensions 20 x 22 x 22. 
In doing this we follow two points: (i) the verification of 
Eq. (|43p and (ii) the influence of the mass increase on 
the observables, especially the shear rate. The geometry 
factor Cg for the simulations can be calculated using Eq. 
(im. The number of time steps in all the simulations is 
Nt — 5 - 10*. Both the inlet / outlet boundary conditions 
proposed by Skordos and by Latt have been tested sep- 
arately. Through the relaxation parameter r the Mach 
number is adjusted. 

The results for the time evolution of the mean density 
and the total mass increase after 5 • 10* time steps as a 
function of the Mach number are presented in figure 21 
As a matter of fact, the time evolution of the mean den- 
sity is well described by a simple exponential (top plot 
in figure 2]) . The recovered exponent for the presented 
example and the theoretical value only differ by 8%. An 
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Table II. Overview of the locations Topt and the relative errors e(ropt) of the minima for different boundary conditions (BC) 
and Reynolds / Mach numbers. 
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Figure 4. Mean density p (top) in lattice units for Re = 5 
and Ma = 0.15 as a function of time step: The simulation 
data (theoretical value CgMa"^/Re = 1.17- 10"*) is excellently 
reproduced by a simple exponential with CgMa^/Re = 1.26 ■ 
10~*. Mass ratios In/i :— ln(A/5oooo/Afioooo) (bottom) for 
different Re as function of the Mach number: The slopes are 
close to the theoretical value of 3. Note that In itself is 
plotted logarithmically. 



even better agreement can be found for the Mach num- 
ber dependence (bottom plot in figured]). We present a 
logarithmic plot of In /i versus the Mach number. /Lt > 
is defined as the ratio of the total mass at time step 5 • 10''' 
and 1 • 10''. Using the mass at t = 1 • 10"* instead of i = 0, 
transient effects are avoided. The predicted exponent of 
3 is nearly obtained in all simulations. As a result, the 
simple model presented in section IIIII is able to describe 
the mass increase during the simulation very well. 

Some remarks should be made at this point. In figure 
m only the results using Latt's scheme are shown. The 
reason is that the results from Skordos' scheme basically 
are equivalent. We have observed, however, that the fi- 
nite difference scheme by Skordos becomes unstable for 
values of r larger than k, 1.8. This is not very alarming, 
since usually LB simulations are performed at t < 1. 
The large values of r have been necessary to increase the 
Mach number for fixed Re and lattice size. Especially 
it has been shown that Skordos' scheme is superior to 
other boundary conditions when simulating flows with 
very large Reynolds numbers [1^. We would also like to 
underline that this analysis only makes sense as long as 
all parameters are in a range, where the LBM and the 
boundary conditions are reliable. Increasing r even fur- 
ther leads to stronger deviations from the Ma'^-law (data 
not shown), but this behavior is of no practical interest. 

The other important observation is that the relative er- 
rors Eq. ([5^ of the velocity and the shear rate do not de- 
pend on time, once steady state has been reached. Here, 
'steady state' is understood in the sense that, despite the 
increase of fluid density with time, the velocity fleld be- 
comes time-independent. The possibility to deflne steady 
state alone shows that the mass increase does not affect 
the velocity. In fact, the errors of the velocity and the 
relevant shear component Sxy remain constant within at 
least flve signiflcant digits. This observation even holds 
if the mass increase becomes huge, as for example in the 
case of Re = 5 and Ma — 0.20, where the mean den- 
sity is p > 5.4 • 10* after 5 • 10* time steps. The mass 
increase is, therefore, effectively harmless in the case of 
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Poiseuille flow, which is a very encouraging result. It is 
an open question whether the mass increase affects the 
velocity and shear rate profiles in time-dependent simula- 
tions, since the distribution of the constantly instreaming 
mass may interfere with the time evolution of the physi- 
cal system. 

V. CONCLUSIONS 

The shear stress — beside the velocity — has become 
an important observable of computer simulations, but an 
analysis of its convergence behavior has still been miss- 
ing. We have shown that the shear stress obtained by 
the lattice Boltzmann method is second order accurate 
and thus behaves similar to the velocity. For a body 
force driven flow, a slope parameter of —1.9 has been 
found. Velocity boundary conditions, on the other hand, 
can spoil this convergence rate, but the effective scaling 
remains clearly better than first order (—1.4 and —1.6 
for Skordos' and Latt's boundary conditions [1^, H^, re- 
spectively). For this reason, it is not justified to denote 
the shear stress a first order quantity. It must be noted 
that the high symmetry of the Poiseuille flow may addi- 
tionally decrease the error of the shear rate, but — from 
a theoretical point of view — the second order nature 
should not be affected by a modified geometry or time- 
dependent flows. This statement remains to be tested. 

Like all numerical computations, lattice Boltzmann 
simulations require the setup of some basic dimensionless 
parameters, including the relaxation time t, the Mach 
number Ma and the numerical grid size N. The choice 
of those parameters has a very significant impact on the 
accuracy of the simulations. Supporting Holdych et al. 
(2^ , we come to the conclusion that it is in general recom- 
mended to use a relaxation parameter r « 0.9. This value 



leads to small velocity and shear errors. The accuracy of 
the simulations also depends on the Mach number, if ve- 
locity boundary conditions are used. Those boundary 
conditions become less reliable at large Mach numbers 
and may spoil especially the results for the shear stress. 
If large Mach numbers cannot be avoided, using a finite 
difference scheme for the shear rate is a good alternative, 
since it is more robust under those circumstances. 

The use of velocity boundary conditions in the lattice 
Boltzmann method can lead to an artificial increase of 
the mass in the numerical grid. For three dimensional 
Poiseuille flow, we have demonstrated theoretically that 
the mass increase is only a function of the Mach number. 
If velocity boundary conditions at both inlet and out- 
let are used, the mass increase can only be reduced by 
choosing a smaller value for the Mach number. From a 
simple theoretical consideration it follows that the mass 
increase scales with Ma^. This result has been confirmed 
with success in our simulations. The increase of the mass 
of the numerical grid does not influence the values of the 
velocity and the shear rate, since their definitions do not 
contain the local density. Although the mass increase 
is not physical, it basically poses no threat to the qual- 
ity of the simulations. However, the velocity boundary 
conditions proposed by Skordos have been observed to 
become numerically unstable, if the relaxation parame- 
ter T is larger than « 1.8. 
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